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We study a genetic regulatory network model developed to demonstrate that genetic robustness 
can evolve through stabilizing selection for optimal phenotypes. We report preliminary results on 
whether such selection could result in a reorganization of the state space of the system. For the 
chosen parameters, the evolution moves the system slightly toward the more ordered part of the 
^— K. ■ phase diagram. We also find that strong memory effects cause the Derrida annealed approximation 

to give erroneous predictions about the model's phase diagram. 

I. INTRODUCTION 

Gene networks are extremely robust against genetic perturbations 0,0]- For example, systematic gene knock-out 
studies on yeast showed that almost 40% of genes on chromosome V have no detectable effects on indicators like cell 
division rate ||, Similar studies on other organisms agree with these results 0,0]. ^ ls a l so known that phenotypically, 
most species do not vary much, although they experience a wide range of environmental and genetic perturbations. 
This striking resilience makes one wonder about the origins, evolutionary consequences, and mechanistic causes of 
* ^ . genetic robustness. 

It has been proposed that genetic robustness evolved through stabilizing selection for a phenotypic optimum. Wagner 
O" 1 showed that this in fact could be true by modeling a developmental process within an evolutionary scenario, in which 
the genetic interaction sequence represents organismal development, and the equilibrium configuration of the gene 
| network represents the phenotype Q. His results show that the genetic robustness of a population of model genetic 
■ regulatory networks can increase through stabilizing selection for a particular equilibrium configuration (phenotype) 
of each network. 

In this paper we investigate the effects of the biological evolution of genetic robustness on the dynamics of gene 
regulatory networks in general. In particular, we want to answer the question whether the evolution process moves 
the system to a different point in the phase diagram. Below, we present some preliminary results. 
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II. MODEL 
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We use a model by Wagner Q, which has also been used by other researchers with minor modifications. Each 
q-( individual is represented by a regulatory gene network consisting of TV genes. The expression level of each gene, s,-, 
has only two values, +1 or — 1, expressed or not, respectively. The expression states change in time according to 
regulatory interactions between the genes. The time evolution of the system configuration represents an (organismal) 
developmental pathway. The discrete-time dynamics are given by a set of nonlinear difference equations representing 
$_i ' a random threshold network (RTN), 



where sgn is the sign function and Wij is the strength of the influence of gene j on gene i. Nonzero elements of 
the N x N matrix W are independent random numbers drawn from a standard normal distribution. (The diagonal 
elements of W are allowed to be nonzero, corresponding to self-regulation.) The (mean) number of nonzero elements 
in W is controlled by the connectivity density, c, which is the probability that a Wij is nonzero. 
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The dynamics given by Eq. ([I]) can have a wide variety of features. For a specified initial configuration s(0), 
the system reaches either a fixed-point attractor or a limit cycle after a transient period. The lengths of transients, 
number of attractors, distribution of attractor lengths, etc. can differ from system to system, depending on whether 
the dynamics are ordered, chaotic, or critical. The fitness of an individual is defined by whether it can reach a 
developmental equilibrium, a certain fixed gene-expression pattern, s*, in a "reasonable" transient time. Further 
details of the model are explained in the next section. 

III. MONTE CARLO SIMULATIONS 
A. Generation and Robustness Assessment of Random Networks 

We studied populations of N = 400 random networks (founding individuals) with N — 10. Each network was 
assigned a matrix W and an initial configuration s(0). W was generated as follows. Each Wij was independently 
chosen to be nonzero with probability c. If so, it was assigned a random number drawn from a standard gaussian 
distribution, N(n = 0, a = 1). Then, each "gene" of the initial configuration, s.;(0), was assigned either —1 or +1 at 
random, each with probability 1/2. 

After W and s(0) were created, the dynamics were started and the network's stability was evaluated. If the system 
reached a fixed point, s*, in 3N timesteps, then it was considered stable and kept. Otherwise it was considered 
unstable, both W and s(0) were discarded, and the process was started over and repeated until a stable network 
was generated. For each stable network, its fixed point, s*, was regarded as the "optimal" gene-expression state 
(phenotype) of the system. This is the only modification we made to Wagner's model: he generated networks with 
preassigned s(0) and s*, whereas we accept any s* as long as it can be reached within 37V timesteps from s(0). 

After generating N = 400 individual stable networks, we analyzed their state-space structures and evaluated their 
robustness as discussed in subsection IIII CI 

B. Evolution 

In order to generate a breed of more robust networks, a mutation-selection process was simulated for all of the 
W = 400 random, stable networks as follows. First, a clan of N 1 — 500 identical copies of each network was generated. 
For each clan, a four-step process was performed for T — 400 generations: 

1. Recombination: Each pair of the N rows of consecutive matrices in the clan were swapped with probability 1/2. 
Since the networks were already shuffled in step 4 (see below), there was no need to pick random pairs. 

2. Mutation: Each nonzero Wij was replaced with probability l/(oZV 2 ) by a new random number drawn from the 
same standard gaussian distribution. Thus, on average, one matrix element was changed per matrix per Monte 
Carlo step. 

3. Fitness evaluation: Each network was run starting from the original initial condition, s(0). If the network 
reached a fixed point, s^, within 37V timesteps, then its fitness, f(st,s*) = exp(— H 2 (s^, s*)/cr s )), was calculated. 
Here H(s^,s*), denotes the normalized Hamming distance between and s* , and a s denotes the strength of 
selection, s* is the optimal gene-expression state, which is the final gene-expression state of the original network 
that "founded" the clan. We used a s = 0.1. If the network could not reach a fixed point, then it was assigned 
the minimum nonzero fitness value, exp(— l/cr s ). 

4. Selection/ Asexual Reproduction: The fitness of each network was normalized to the fitness value of the most 
fit network in the clan. Then, a network was chosen at random and duplicated into the descendant clan with 
probability equal to its normalized fitness. This process was repeated until the size of the descendant clan 
reached N' . Then the old clan was discarded, and the descendant clan was kept as the next generation. Note 
that this process allows multiple copies (offspring) of the same network to appear in the descendant clan, while 
some networks may not make it to the next generation due to genetic drift. 

At the end of the T — 400 generation selection, any unstable networks were removed from the evolved clan. 
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FIG. 1: (a) x{t + 1) shown vs. x(t) for N = 16 and (fc) = 4. The theory, Eq. ©, is in good agreement with the simulations. 
The deviations are due to the small size of the simulated system as the theoretical calculation assumes N ^> (k). (b) Damage- 
spreading rate, x(t + 1) — x(t) vs. x{t), for random and evolved networks with N — 10 and (fe) = 5 and 7, showing the 
difference between the "random" and "evolved" curves. Only the first half of the curves are shown since x(t + 1) vs. x(t) is 
point-symmetric about (1/2, 1/2). The results were averaged over 10 random networks and all of their evolved descendants 
(~ 300 evolved networks per random network). The evolved curves for each (k) lie very close to their "random" counterparts. 
However, they are outside twice the error bar range of each other at most data points. 



C. Assessment of Robustness 

The mutational robustness of a network was assessed slightly differently for random and evolved networks. For a 
random network, first, one nonzero toy was picked at random and replaced by a new random number with the same 
standard gaussian distribution. Then, the dynamics were started, and it was checked if the system reached the same 
equilibrium state, s*, within 3N timesteps. This process was repeated 5000c times using the original matrix (i.e., 
each mutated matrix was discarded after its stability was evaluated) . The robustness of the original network before 
evolution was defined as the fraction of singly-mutated networks that reached s*. 

For the evolved networks, clan averages were used. For each of J\f opt < 400 networks in a clan, robustness was 
assessed as described above with one difference: the number of perturbations was reduced to 5000c/A/ r ° pt per network 
to keep the total number of perturbations used to estimate robustness of networks before and after evolution approx- 
imately equal. The mean robustness of the those A/" opt networks was taken as the robustness of the founder network 
after evolution. Therefore, the robustness of a network after evolution is the mean robustness of its descendant clan 
of stable networks. 



IV. RESULTS 

As Wagner pointed out, the stabilizing selection described above increases the robustness of the model population 
against mutations However, it is not very clear what kind of a reorganization in the state space occurs during the 
evolution. Also, it is not known whether this robustness against mutations leads to robustness against environmental 
perturbations. In this paper, we focus on the effects of evolution in terms of moving the system to another point in 
the phase diagram. In other words, we investigate whether the system becomes more chaotic or more ordered after 
evolution. 

A standard method for studying damage spreading in systems such as the one considered here is the Derrida 
annealed approximation 0, @|, in which one calculates changes with time of the overlap of two distinct states, s(t) 
and s(t), 

1 N 

*(<) = 2JV !>*(*) +*(*)!• (2) 

i=l 
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The change of the overlap over one time step for N 3> (k) = Nc is given by 



x(t + 1) = n(0)x(t) 



n(l)x(t) +2_^n(k) 

k=2 



k-1 



(x(t)) k +J2^k(l)V(k,l) 



1=1 



(3) 



where the Poisson distribution n(k) = n k cxp(— n)/kl, is the probability of finding a gene, i, with k inputs, the binomial 
distribution Tlk{l) = ( ; ) (1 — x(t)) 1 (x(t)) k ~ l is the probability of finding k — l of these inputs in the overlapping parts 



of s(t) or s(t), and V(k,l) = 1 — -|arctan \^\/l/{k — I) J (for k > I) is the probability of the sum of k — l matrix elements 

being larger than the sum of I matrix elements, which are independent and iV(0, 1) distributed. Here, k = (k), the 
mean number of inputs per node. 

For most RTNs that have been studied so far [1,1, Eq. ® can be iterated as a map to give the full time evolution of 
the overlap. Changes in the fixed-point structure of this map with changing (k) would then signify phase transitions 
of the system. As seen in Fig. QJi, for (k) = 4, such a map would have a stable fixed point at x = 1/2. One 
can also show that lim x ^^ + dx(t + l)/dx(t) > 1 for all (k) > (this implies lim,^)^!- dx(t + l)/dx(t) > 1 and 
lim x (t)^i/2 dx(t+l)/dx(t) < 1), and so it would seem that the system has no phase transition and always stays chaotic 
for nonzero (k). However, simulations of damage spreading for longer times [7] indicate that the system studied here 
has strong memory effects due to the update rule for spins with no inputs, given by the last line in Eq. ([1]), which 
retard the damage spreading [§]. In fact, like other RTNs the system undergoes a phase transition near (k) « 2 
from a chaotic phase at larger (k) to an ordered phase at smaller (fc). The strong, retarding memory effects mean 
that Eq. ([3]) cannot be iterated as a map, and the naive prediction based on the Derrida annealed approximation is 
erroneous. 

Despite its irrelevance for the long-time damage spreading, the damage-spreading rate shown in Fig. [TJd properly 
describes the short-time dynamical character of the system. However, as Eq. assumes that the interaction constants, 
Wij , are statistically independent, it may not apply to evolved networks as we do not know whether the selection process 
creates correlations between the matrix elements. Nevertheless, we can still compute x(t + 1) as a function of x(t) 
numerically to see if there is a change in the degree of chaoticity (or order) of the dynamics. As seen in Fig. [TJd, the 
damage-spreading rates for evolved networks are slightly (but statistically significantly) lower than for their random 
predecessors, which are thus slightly more chaotic. 

To summarize, we have presented preliminary results on some general properties of a popular RTN model of a gene 
regulatory network and on how the biological evolution of genetic robustness affects its dynamics [§]. We have also 
shown that the update rule for spins without inputs leads to strong memory effects that invalidate naive iteration 
of the Derrida annealed approximation as a map. The evolutionary process that improves the genetic robustness of 
such networks has only a very small effect on their dynamical properties: after evolution, the system moves slightly 
toward the more ordered part of the phase diagram. 
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